function xdot = x_dot(sc, u_opt, auxdata)

mu      = auxdata.mu;
Tmax    = auxdata.Tmax;
c       = auxdata.c;

xdot = zeros(7,1);

p = sc(1);
f = sc(2);
g = sc(3);
h = sc(4);
k = sc(5);
L = sc(6);

q   = 1 + f*cos(L) + g*sin(L);
s2  = 1 + h^2 + k^2;               % angular momentum magnitude

A(1,1) = 2*p/q * sqrt(p/mu);
A(1,2) = 0;
A(1,3) = 0;
A(2,1) = sqrt(p/mu)/q * ( (q+1)*cos(L) + f );
A(2,2) = sqrt(p/mu) * sin(L);
A(2,3) = - sqrt(p/mu) * g/q * ( h*sin(L) - k*cos(L) );
A(3,1) = sqrt(p/mu)/q * ( (q+1)*sin(L) + g );
A(3,2) = - sqrt(p/mu) * cos(L);
A(3,3) = sqrt(p/mu) * f/q * ( h*sin(L) - k*cos(L) );
A(4,1) = 0;
A(4,2) = 0;
A(4,3) = sqrt(p/mu) * s2 * cos(L) / (2*q);
A(5,1) = 0;
A(5,2) = 0;
A(5,3) = sqrt(p/mu) * s2 * sin(L) / (2*q);
A(6,1) = 0;
A(6,2) = 0;
A(6,3) = sqrt(p/mu) * ( h*sin(L) - k*cos(L) ) / q;

%% perturbation
fp = perturbation(sc, auxdata);

%% control modified
xdot(1:6,1) = A * (u_opt + fp) + [0;0;0;0;0;sqrt(mu*p)*(q/p)^2];
xdot(7,1)   = - Tmax/c;
xdot = real(xdot);


















